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ABSTRACT 


This  memorandum  examines  an  application  of  the  least-squares  approximation  of  certain 
structured  matrixes  to  a  new  algorithm  for  the  maximum  likelihood  (ML)  estimation  of 
directions  of  arrival  (DOAs)  and  frequencies  of  sinusoids.  The  method  involves  reordering  the 
steps  in  the  standard  ML  search  algorithm  so  that  the  majority  of  the  calculations  are  made  prior 
to  execution  of  the  multidimensional  search.  The  computational  efficiency  is  due  to  the  fact  that 
only  scalar  operations  are  performed  during  the  search.  This  new  fast  maximum  likelihood 
estimation  (FMLE)  algorithm  is  more  than  one  order  of  magnitude  faster  than  the  standard 
direct-search  ML  method. 
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A  Fast  Maximum  Likelihood  Multidimensional 

Search  Algorithm 


1.  INTRODUCTION 


The  maximum  likelihood  (ML)  estimation  method  is  not  always  the  method  of  choice 
for  many  applications  because  of  the  complexity  of  the  calculations  involved.  This  memo¬ 
randum  describes  a  more  efficient  computational  method  that  will  allow  ML  estimation  to 
be  used  for  a  wider  class  of  problems  than  at  present. 

The  fast  maximum  likelihood  estimation  (FMLE)  algorithm  introduced  here  can  be 
applied  to  two  problems:  direction-of-arrival  (DOA)  estimation  and  sinusoidal  frequency 
estimation.  Under  mild  restrictions  on  the  statistics  of  the  noise,  these  problems  are  shown 
to  be  equivalent  to  a  structured  least-squares  problem.  An  estimate  of  a  noise-corrupted  ma¬ 
trix  of  data  can  be  improved  by  taking  advantage  of  prior  knowledge  concerning  its  under¬ 
lying  structure  (for  example,  the  underlying  Toeplitz  structure  of  a  covariance  matrix).1,2,3 
Another  type  of  estimate  is  one  that  is  both  structured  and  low  rank.4’5 

In  this  study,  the  ML  matrix  estimate  is  required  to  have  a  specific  rank,  which  is  equal 
to  the  number,  N,  of  sources  present.  In  addition,  the  estimate  is  required  to  have  a  column 
space  spanned  by  a  set  of  N  basis  vectors  that  have  a  specific  structure.  In  this  manner, 
the  ML  problem  may  be  thought  of  as  a  type  of  structured  low-rank  approximation. 

In  section  2,  two  structured  least-squares  approximation  problems  are  introduced  and 
the  solutions  are  written  in  a  form  involving  the  maximization  of  a  single  expression. 
Section  3  examines  applications  of  ML  estimation  to  the  determination  of  DOAs  and  to 
the  estimation  of  frequencies  of  sinusoids  in  noise.  In  section  4,  the  FMLE  method  is 
developed.  After  the  results  of  the  simulations  axe  given  in  section  5,  section  6  concludes 
the  memorandum. 


2.  FORMULATION  OF  TWO  LEAST-SQUARES  PROBLEMS 


In  this  memorandum,  the  columns  {dl5  •  •  • ,  dr}  of  the  matrix  D,  over  which  each  search 
is  to  be  made,  are  of  a  special  form.  Each  column  of  D  can  be  written  as  an  L- vector: 

t?  =  [l  z  z2  -  zi-1]T. 

The  vector,  v,  will  be  called  a  Fourier  vector,  because  when  2  =  etS ,  the  inner  product  vHw 
can  be  interpreted  as  the  discrete  Fourier  transform  (DFT)  of  w,  evaluated  at  2.  For  DOA 
problems,  Fourier  vectors  are  also  referred  to  as  steering  vectors.  The  first  least-squares 
problem  is  presented  next: 


1 


Problem  1—  Let  A  be  a  complex  L  x  s  matrix.  Find  an  L  x  r  matrix ,  D  =  [</1?  ■  •  ■ ,  dr\, 
with  distinct  Fourier  vectors  as  columns  and  an  r  x  s  matrix,  p,  such  that  the  following  is 
minimized: 

argmin  \\A- Dp\\F  .  (1) 

U,p 

We  then  have 

arg  min  \\A  —  Dp\\2F  =  arg  min  trace  (.AH A  —  AH Dp  —  pH Du  A  +  pH DH Dp)  .  (2) 

D,p  D,p  '  J 

Using  the  projection  theorem,6  the  columns  of  the  error  matrix  A  —  Dp  of  the  ap¬ 
proximation  are  orthogonal  to  the  column  space  in  which  the  solution  lies.  Therefore,  the 
orthogonal  projection  of  the  error  onto  the  column  space  of  D  will  be  zero: 


D(DhD)~1Dh  (A  -  Dp)  =  0  . 

Since  D  is  a  Vandermonde  matrix  with  distinct  columns,  DH D  is  nonsingular.7  We  then 
have 

p  =  (DHD)~lDHA  . 

From  equation  (2),  it  follows  that 

argmin  \\A  —  Dp\\F  —  arg  max  trace  (A11 D{Dh  D)~l  DH  A')  .  (3) 


Problem  2  —  Let  Abe  a  complex  Lx  L  Hermitian  matrix.  Find  an  Lxr  matrix,  D,  where 
the  columns,  di,i  —  1,  -  •  • ,  r,  of  D  are  Fourier  vectors  and  an  r  xr  matrix,  p,  such  that  the 
following  is  minimized: 


arg  min 

D,p 


A  -  DpDH 


(4) 


In  a  similar  manner  to  the  previous  case,  this  leads  to  the  same  minimization  problem: 
argmin  \\A  —  DpDH\\2F  =  arg  max  trace  ( AHD(DHD)~1DHA )  ,  (5) 

D,p  O'  / 

where  p  =  (DH D)-1  DH  AD(DH D)~l .  Let  D  be  any  orthonormal  set  of  vectors  spanning 
the  same  column  space  as  D.  Due  to  the  uniqueness  of  an  orthogonal  projection  matrix, 
we  have 

D(DhD)~1Dh  =  DDh. 

We  can  then  further  simplify  the  right-hand  sides  of  equations  (3)  and  (5)  to 

aigmax\\DH A\\2F  .  (6) 

The  least-squares  estimates  considered  above  and  the  ML  estimation  methods  to  be 
considered  in  the  next  section  are  both  related  to  a  third  type  of  estimate,  which  is  the 
minimum  variance  unbiased  estimate.  If  the  noise  is  Gaussian  and  the  error  covariance  is 
the  identity  matrix  scaled  by  the  noise  variance,  then  these  three  estimates  axe  equivalent.6,8 
We  assume  that  these  restrictions  hold. 
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3.  DEVELOPMENT  OF  THE  MAXIMUM  LIKELIHOOD 
ESTIMATION  EQUATIONS 


3.1  ESTIMATION  OF  DIRECTIONS  OF  ARRIVAL 


This  section  examines  an  application  of  ML  estimation  to  the  determination  of  DOAs 
in  an  array  processing  scenario,  where  the  signals  may  be  temporally  coherent.  This  case 
cannot  be  dealt  with  directly  by  subspace-based  methods  without  further  processing  9  or 
array  partitioning.10 

The  DOA  problem  treated  here  arises  from  M  narrowband  farfield  signals  impinging 
on  a  linear  array  of  L  equally  spaced  sensors.  The  sensor  outputs  are  subject  to  additive 
zero-mean  white  Gaussian  noise.  The  equation  describing  this  scenario  is 

x(k)  =  B(k)s(k)  +  n(lfc),  k=  1,2, ,  (7) 

where  x(Ar),  called  a  data  snapshot,  is  an  L-vector  of  array  outputs  at  sampling  time  k. 
The  columns  of  B(k)  =  [6i(&)  fh(k)  •  •  •  are  the  steering  vectors  at  time  k: 

bi(k)  =  [1  ex.j>(2nj(3  cos(0t(fc)))  •  •  •  exp(27rj(£  —  1  )/3  cos(0,-(fc)))]r  ,  i  =  1,  •  •  • ,  M. 

Here,  0,-  is  the  ith  bearing  relative  to  broadside  and  (3  is  the  interelement  spacing  of  the 
array  in  wavelengths.  The  value  of  M  can  be  estimated,  for  example,  by  using  the  model- 
based  approach  given  in  reference  11.  The  elements  of  the  signal  vector,  s(k),  are  complex 
exponentials: 

3i(k)  =  a(k)  exp (j<pi(k)),  i  =  1,  ■  •  ■ ,  M,  k  =  1, 2,  •  •  •  . 

Let  E{x}  be  the  expected  value  of  the  random  vector,  x.  We  then  have 

C  =f  E  {xx*}  =  E  {(Bs  +  n)(Bs  +  n)*}  =  B  E  {ss*}  BH  +  E  {nn*}  =  BSBH  +  a2 1  . 

The  ML  estimate  C  of  the  estimated  covariance  matrix  is 

C  =  Tf  £*(*)*(*)*  •  (8) 

i=l 

It  is  assumed  here  that  the  source  DOAs  are  constant  during  the  time  period  over  which 
the  average  in  equation  (8)  is  taken.  The  solution  of  equation  (7)  corresponding  to  the 
concentrated  likelihood  function  in  reference  8  is 

arg  max  (* (»)"(/>(£>H.£>)-‘0,,)z(i))  ,  (9) 
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where  x(i)  is  the  data  snapshot  obtained  at  time  i.  The  ML  estimation  of  the  true  covariance 
matrix,  C,  can  be  obtained  by  taking  N  data  snapshots  and  averaging  the  expression  in 
equation  (9)  over  the  following  terms: 

1  N 

—  arg  max  trace  (x(i)H  D(DH  D)~l  DH  x(i)} 

=  arg  max  trace  ^D(DH D)~l  PH  ^ 

=  arg  max  trace  (jD(DH  D)~l  DH C'j 
=  arg  max  trace  (UH D(DH D)~x DH U)  , 

where  UUH  is  the  Cholesky  factorization  of  C. 

3.2  ESTIMATION  OF  FREQUENCIES  OF  SINUSOIDS 


Consider  taking  L  consecutive  samples  from  a  time  series  containing  a  weighted  sum  of 
M  exponentials  in  Gaussian  white  noise: 

M 

y(n )  =  a*'  exp(27rj/,n)  +  w(n),  n  =  1,  •  ■  • ,  L  . 

i= 1 


The  ML  method  applied  to  this  problem12  involves  the  minimization  of  the  following 


error: 


L 

e  =  E 

71=1 


M 

y(n)  -  exP  fajfcn) 


i=i 


where  the  a^s  and  f[s  axe  to  be  estimated.  Let  Y  =  [y (1)  •  •  •  y(L)]T.  The  least-squaxes 
solution  then  is 

aigmax(YHD(DHD)~1DHY)  ,  (10) 

from  which  a  squaxe  root  can  be  calculated:  (DH D)~l  =  F2.  It  follows  that  the  columns 
of  DF  axe  orthonormal.  With  this  result,  equation  (10)  can  be  rewritten  as 


arg  min  \\FH DHY\\p  . 


(11) 


An  alternate  orthonormalization  for  D  is  considered  in  the  next  section. 


4.  DERIVATION  OF  THE  FMLE  ALGORITHM 


We  now  describe  a  way  to  take  advantage  of  the  simplified  form  of  the  solution  in  the 
expression 

axgm^.x||T>'ffA||^  .  (12) 
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In  the  following  discussion,  assume  that  A  is  a  given  LxQ  matrix  and  that  D  =  [d\  -  -  -  dr\ 
is  an  L  x  r  matrix  (r  <=  Q),  with  distinct  Fourier  vectors  as  columns. 

The  columns  of  D  are  to  be  selected  from  a  given  L  x  K  matrix  B  of  Fourier  vectors, 
and  the  multidimensional  search  is  performed  over  all  subsets  of  r  distinct  columns  of  B. 
We  seek  an  alternative  to  using  the  factor,  F,  in  equation  (11)  for  orthonormalizing  the 
columns  of  D.  Since  the  Frobenius  norm  is  unitarily  invariant,  all  orthonormal  bases  for  the 
column  space  of  D  are  equivalent,  with  respect  to  the  value  of  the  expression  in  equation 
(12).  We  will  use  the  Gram-Schmidt  procedure  to  accomplish  the  orthonormalization  of 
the  columns  of  D. 

In  order  to  calculate  the  inner  products  dfdj  between  the  Fourier  vectors  d,  and  dj,  we 
define  the  matrix  [7]: 

It  will  be  assumed  that  the  columns  of  B  are  normalized;  ( BHB)a  =  1,  i  =  1,  •  •  • ,  K.  It  is 
easily  shown  that  BH B  is  Toeplitz,  as  well  as  Hermitian.  Therefore,  only  the  first  row  of 
BH B  needs  to  be  formed.  Consider  two  Fourier  vectors,  dn  and  dTO,  where  dn  is 

dn  =f  [1  exp (jfin)  •  •  -  exp {jn(L  -  l)n)]r  . 

We  now  have 

„  Wjff,  1-exp  (-jft(n-m)L) 

n  m  j  _  exp (— j/i(n  -  m)) 

where  j  =f  y/—l.  We  also  define  the  matrix  [a]  =f  BH A,  where 

Onro  =  (fi  •^•)nTF»5  ^  =  l)‘"i  K ,  W  =  1,"  ■  •  ,  Q  . 

The  inner  products  in  this  matrix-matrix  product  involve  the  unorthogonal  Fourier  vectors. 
We  next  describe  the  search  algorithm.  To  illustrate  the  overall  procedure,  let  r  =  2.  We 
then  execute  the  two  nested  loops  as  follows:  for  i  —  1,  •  •  ■ ,  M  —  1;  for  j  =  i  +  1,  •  •  • ,  M. 

Consider  the  entries  of  the  matrix  DHA1  where  D  —  [di  dj].  The  maximization  is 

argmax  \\D*H A\\2F  , 

where  D*  is  any  orthonormalization  of  the  columns  of  D.  Given  i  and  y,  the  sum  of  the 
squared  magnitudes  of  the  terms  in  the  first  row  of  DH  A  is 

£  \d?ak\2  , 

fc=i 

where  a k  is  the  kth  column  of  A.  For  the  second  row,  we  must  use  the  orthonormalized 
values  for  dj  given  di: 


(14) 


(13) 


4  ®  (4  -  MfaKVAi,  i  =  i.' 
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where 


Ai  =  IKil|2  =  |K-(^K-||2 

is  the  normalizing  factor.  The  norm  j3{j  of  d3  can  be  found  from  the  matrix  [7]  since 


0 if  -  0 ~  (dpi)*) 

1  -  |7«|2 


(15) 


for  i  ^  j.  In  the  general  case  for  r  >  2,  equation  (15)  becomes 

ip,  =  U-Y,{dfd,)d^  U-f.iip,)^ 

=  1  —  2  ^  |7ij  |  +  ^  7nj7mj7nm  • 

i=l  n=l  m=l 

The  sum  of  the  squared  magnitudes  of  the  elements  in  the  second  row  for  D^H  A,  given 
di ,  is  then 

E  •  (16) 

fc=  1 

Consider  the  following  term  from  equation  (16): 

dfak  =  ((d?-(d?dj)d?)/0ij)ak 
=  (dfak  -  Jijdfak)  /fa 

=  (ay*  -  7*yQ*>)  /Ay,  *  =  !»•••»(?•  (17) 

In  this  key  equation,  the  computational  savings  to  be  realized  are  due  to  the  fact  that  all 
the  individual  terms  are  scalars,  and  each  term  can  be  calculated  before  the  execution  of 
the  nested  loops  constituting  the  search.  Combining  equations  (14),  (16)  and  (17),  the 
expression  for  ||Z}t'ff/4||jr,  where  D  =  [di  dj],  is 


\ 


Wa^f 


H  ' °*l2 


fc=l  *=1 


lQ»fc|2  +  Y1  \art  -  7tya»*l2  /Pij  • 


*=1  Jt=l 


(18) 


It  can  be  seen  from  the  simulations  that  this  approach  results  in  over  an  order  of 
magnitude  gain  in  efficiency  when  compared  with  the  standard  ML  search  procedure.  In 
summation  of  this  section,  the  newly  developed  FMLE  algorithm  for  an  r-dimensional 
search  procedure  involves  the  following  steps: 

1.  Compute  the  distinct  entries  of  BH B  at  the  beginning  of  the  analysis,  using  equation 
(13). 
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2.  Compute  BH A  before  entering  the  nested  loops.  This  can  be  done  efficiently  with  an 
update  procedure  whenever  a  sliding  window  or  an  exponential  window  is  used.  If  the 
Fourier  vectors  in  the  columns  of  B  are  equally  spaced  in  angle,  BH  A  can  be  found 
by  computing  the  DFTs  of  the  columns  of  A  and  summing  the  squared  magnitudes 
of  the  elements  in  each  row. 

3.  Conduct  the  r-level  multidimensional  search. 

The  computational  complexity  of  FMLE  is  0(Kr),  where  there  is  an  r-level  multidi¬ 
mensional  search  over  a  set  of  K  Fourier  vectors.  The  standard  ML  search  for  this  problem 
is  also  of  0(Kr)  complexity,  but  the  constant  is  much  larger. 

5.  SIMULATIONS 


In  the  first  simulation,  as  shown  in  figure  I,  the  performance  of  the  ML  and  FMLE 
algorithms  is  compared  with  the  number  of  array  elements.  Since  there  are  only  scalar 
operations  performed  within  the  FMLE  search,  the  FMLE  complexity  is  of  lower  order 
than  the  ML  complexity. 


Figure  1:  Comparison  of  Standard  ML  and  FMLE  Performance  Versus  Array 
Size. 


The  second  simulation,  as  shown  in  figure  2,  compares  the  ML  and  FMLE  performance 
versus  the  number  of  sources  present,  i.e.,  the  dimension  of  the  multidimensional  search. 


Figure  2:  Complexity  Versus  Dimension  of  Search. 


6.  CONCLUSIONS 


It  has  been  shown  that  the  ML  estimation  of  DOAs  or  frequencies  of  sinusoids  can  be 
formulated  as  a  least-squares  matrix  approximation  problem,  under  mild  restrictions  on 
the  noise  structure.  The  efficiency  of  the  FMLE  algorithm  is  made  possible  by  the  perfor¬ 
mance  of  all  nonscalar  calculations  before  the  execution  of  the  nested  loops  constituting 
the  multidimensional  search.  This  approach  results  in  an  order  of  magnitude  speed-up  in 
the  ML  search. 
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